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Large-eddy simulation of a plane wake 


By S. Ghosal AND M. Rogers 1 


1. Motivation and objectives 

In a previous report (Ghosal et a/., 1992, 1994) the theoretical development 
leading to the dynamic localization model (DLM) for large-eddy simulation (LES) 
was presented. The method has been successfully applied to isotropic turbulence 
(Ghosal et a/., 1992, 1993, 1994, Carati et a/., 1994 - see also the report in this 
volume), channel flow (Cabot, 1993 - see also the report in this volume) and the 
flow over a backward-facing step (Akselvoll & Moin, 1993a & b). Here we apply 
the model to the computation of the temporally developing plane wake. The two 
main objectives of this project are: 

(A) Use the model to perform an LES of a time developing plane wake and com- 
pare the results with direct numerical simulation (DNS) data to see if important 
statistical measures can be reliably predicted. Also, to provide a relative evaluation 
of the several versions of the model in terms of predictive capability and cost. 

(B) If the tests in (A) show that the model generates reliable predictions, then 
use the LES to study various aspects of the physics of turbulent wakes and mixing 
layers. 

According to the notation introduced earlier (see the references above), we rec- 
ognize four versions of DLM: 

(1) Dynamic model (DM): Special case of the DLM applicable only to flows with 
homogeneous directions. 

(2) Dynamic localization model (constrained) [DLM(+)]: A limited version of the 
more general DLM, explicitly prevents backscatter by enforcing a positivity require- 
ment on the Smagorinsky coefficient. 

(3) Dynamic localization model (k-equation) [DLM(k)]: Extended version of DLM 
that incorporates backscatter by introducing a budget equation for the sub-grid 
kinetic energy. 

(4) Dynamic localization model (stochastic) [DLM(S)]: Alternate extension of DLM 
that incorporates backscatter by a stochastic term. 

Tests of the DM and DLM(+) will be presented in this report. The more elaborate 
models DLM(k) and DLM(S) that incorporate backscatter have not yet been tested 
for this flow. In the next section we briefly review the two versions of the model 
tested. No derivations are presented here; the reader is referred to the appropriate 
references (Ghosal et a/., 1992, 1994, and references therein) for the underlying 
theory. 

1 NASA Ames Research Center 
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2. Accomplishments 


2.1 Background 


2.1.1 The dynamic model (DM) 

In its present form, the dynamic model can be written in the following way. For 
homogeneous turbulence the coefficient C is constant in space (but it may be time 
dependent) and is given by 

Here L tJ _ = u t Uj —u t Uj is the Leonard term and m*j = A 2 |5|5,j — A 2 |S|Sjj, where 
Hi and S tJ are the filtered velocity and strain rates and the ‘hat’ denotes the “test 
filtering” operation: 

/(*) = /c(x, y )/( y )dy. (2) 

The ‘grid-level’ filter-width is A (usually taken to be of the order of the grid spac- 
ing) and A (A > A) is the ‘test-level’ filter-width. The angular brackets denote 
averaging over the volume of the domain. 

For flows that are not completely homogeneous but have one or two homogeneous 
direction(s) the DM can still be applied provided one assumes that the u test filtering” 
operation is performed only in the homogeneous direction(s). Such an assumption 
can be justified if the grid in the inhomogeneous direction(s) is so fine that the flow is 
fully resolved in that direction, but in general it is not strictly valid. If one considers 
a flow (such as the plane wake considered in this report) that is homogeneous in 
the x — z plane but inhomogeneous in y, then the DM can be written as 


C(y,t) = piMlL. 


( 3 ) 


where the angular brackets now denote averaging over the homogeneous x — z planes. 

A serious problem with the DM is that it can be applied only to homogeneous 
flows or (under additional assumptions) to flows with at least one homogeneous 
direction. This deficiency is removed by the DLM described next. 

2.1.2 The dynamic localization model: constrained [DLM(-h)] 

In DLM(+) one obtains C(x) as a function of position at each time-step by solving 
an integral equation 


C(x) = 


/(*)+ / £(x>y)C(y)dy 
J J + 

where the suffix indicates the positive part and 

/( x ) = afc;(x)a»(x) q *>( x )-M x ) “ /M x ) J Lij(y)G(y,x)dy 


( 4 ) 


( 5 ) 
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N /CA(x,y) + /CA(y,x)-/C s (x,y) 

/C(X ’ y } " a«(*)a„(x) ’ W ’ 

/C A (x,y) = a ij (x)/3 ij (y)G(x,y) (7) 

and 

£s(x,y) = 0ij(x)0ij( y) J G(z,x)G(z,y) dz. (8) 

In these expressions G(x, y) is the “test filter” a, ; = — 2A 2 |S'|Sj J , flij = — 2A 2 |S , |S; J 
and Lij is the Leonard term. 

The principal weakness of DLM(+) (as well as the DM) is that the restriction 
of C to only positive values is somewhat contrived because it does not account for 
backscatter. However, unlike the DM, the DLM(+) is completely general and can 
be applied to arbitrary inhomogeneous flows. 

2.1.8 The problem of the temporally developing wake 

In a temporally developing wake the flow is statistically homogeneous in the 
streamwise (x) and spanwise ( z ) directions and inhomogeneous in the normal (y) 
direction. The governing equations are the incompressible Navier-Stokes equations 
with periodic boundary conditions in x and z. In the y-direction the domain is 
infinite and the velocity field is assumed to asymptotically approach the free-stream 
velocity, which can be taken as zero in a suitably chosen reference frame. This 
system can be considered to be an approximation to the physically more interesting 
spatially developing wake. If one imagines a 'box’ being advected downstream 
at the Tree-stream’ velocity, then the motion of the fluid in the imaginary box 
approximates a temporally developing wake. The integrated mass flux deficit 





6U(y)dy 


(9) 


is conserved in a temporally developing wake, as opposed to the momentum flux 
deficit 


V* 



+ 6U(y))6U(y)dy , 


( 10 ) 


which is conserved for a spatially developing wake. Clearly, if the mean velocity 
deficit 6U is small compared to the free stream velocity Uoo, then /i* « 17oo/n 
A suitable scale for the velocity is the initial centerplane velocity deficit 6U 0 = 
-(6Z7(0))t=o and a suitable length scale is then n/6U 0 . The corresponding time 
scale is fi/(8Uo) 2 . We will quote most of our results in these units. 


2.2 Computational methods 

The numerical method used is a spectral method in vorticity variables. Both the 
velocity and vorticity are periodic in the x and z directions and can therefore be 
expanded in a basis of trigonometric functions for these variables. The y-direction 
is somewhat more difficult to deal with since the domain is infinite in y. One 
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method is to choose a basis of functions that have an infinite support (such as 
the Jacobi polynomials coupled with a mapping to the infinite interval) for the 
y-direction (Spalart et al., 1991). However, here we use a trick that leads to a 
simpler code. VVe take advantage of the fact that in a wake the vorticity field is 
much more confined in the y direction than the velocity field. One then expands 
the vorticity in a trigonometric series in y defined over (y min, y max) with periodic 
boundary conditions. This is permissible provided that the vorticity is narrowly 
confined around y = 0 and effectively decays to zero at the boundaries y min and 
2/max • The velocity field is not so confined and cannot be represented in terms 
of these trigonometric functions. But once the vorticity field is determined, the 
correct velocity field may be obtained by adding a potential “correction” to the 
periodic velocity field associated with the vorticity field so as to match the boundary 
conditions at y = ±oo. Further details of the computational method may be found 
in Corral and Jimenez (1993). The method of solving the integral equation to 
determine the coefficient C has been described elsewhere (Ghosal et al., 1992, 1994). 
The test filter width in these computations was taken to be twice the grid-filter 
width, A = 2A, and a ‘top-hat’ filter was used with a Simpson’s rule quadrature. 

For initial conditions we take two realizations of ‘turbulence over a flat plate’ 
from DNS data generated by Spalart (1988) and ‘sandwich’ them to produce a 
wake. Physically this corresponds to a situation where two independent boundary 
layers exist on either side of a rigid plate and the plate is instantaneously “dis- 
solved” without disturbing the surrounding fluid. All the parameters in the LES 
are chosen so as to correspond to the “unforced wake” case of Moser and Rogers 
(1994) mentioned above. 

The LES reported here was performed on a grid of size Nx = 64, Ny = 48, and 
N z = 16. Therefore, all DNS data must first be ‘filtered’ to the same resolution 
as the LES. This is done by truncating the DNS data in Fourier space to the 
same number of modes retained in the LES. This filtering procedure is applied to 
the initial conditions as well as to all DNS data with which we wish to compare 
our LES results. The ‘filtered DNS’ represents the theoretical best that can be 
achieved by any LES. The LES with DM took about 11 minutes of CPU time. For 
the DLM(-f) the CPU time depended on the level of convergence required for the 
solution of the integral equation. We measured the degree of convergence by the 
rms error in satisfying the integral equation normalized by the maximum value of 
(C) xz' When it was required that the error as defined above should not exceed 
10“ 4 , the DLM( + ) used about 18 minutes of CPU time. To test if this level of 
convergence was adequate, the simulation was rerun with the convergence criterion 
set at 10 9 . There were no observable differences in any of the computed statistical 
measures. For comparison, the high resolution DNS of Moser and Rogers of the 
same flow over the same time interval cost about 200 CPU hours. All computations 
were performed on a CRAY C90. 


2.S Results 

The gross features of the wake are characterized by the maximum wake deficit 
SU m of the mean velocity profile and the ‘half-width’ b of the wake. The half- width 
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t(SU 0 ) 2 /fi 

FIGURE 1. The square of the wake-width as a function of time using DLM(-i-) 
; DM - - No model ; filtered DNS • 



t(6U 0 ) 2 /n 

FIGURE 2. The product of the wake- width and the maximum velocity deficit as a 
function of time using DLM(+) ! DM — — — i No model ; filtered DNS • 

is defined here as the distance between the two points at which the mean velocity 
deficit is half its maximum value. Fig. 1 shows b 2 plotted as a function of the 
time t for the LES, filtered DNS, and LES with the subgrid model turned off. The 
prediction of the DM is closest to the filtered DNS. The width grows as b ~ \ft 
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y/i> 

Figure 3. The mean wake velocity deficit in self-similar coordinates using 
DLM(+) DM No model ; filtered DNS • 

in the self-similar region (t(6U 0 ) 2 /fi as 50 - 100) as expected. Fig. 2 shows the 
product b(6U m ) as a function of t. In all cases this quantity exhibits a plateau 
during the self-similar period. Note that the Reynolds number Re\, = b 6U m /;/ = 
20006 6U m /fj, as 2000 in the self-similar period. 

Fig. 3 shows the mean velocity profile plotted in self-similar coordinates 6U* - 
6U/6U m and y, = y/b for t (6U 0 ) 2 /fi ss 50 — 100. In all cases very good self-similar 
collapse is observed (even with the subgrid model turned off!). Thus, the mean 
velocity profile is quite insensitive to the subgrid model. 

Figs. 4 (A), (B), (C), and (D) show the second-order velocity statistics (u 2 ), 
(u 2 ), (w 2 ), and < uv ) respectively. Here u, v, and w are the velocities in the x, y, 
and 2 directions, respectively, with the mean velocity subtracted out. The angular 
brackets denote averaging over x-z planes. In all cases it is observed that both the 
DM and DLM(+) predict the second-order statistics very well. The quality of the 
predictions deteriorates significantly if the model is turned off (except for (uv)). The 
better agreement for the (m;) profile is to be expected since it is directly related 
to the mean velocity profile 6U(y) through the x-component of the momentum 
equation and we have already seen that 6U(y) is insensitive to the subgrid model. 

Figs. 5 (A), (B), (C), and (D) show the second-order vorticity statistics (w 2 ), (u> 2 ), 
(uf z ), and (ix> x u>y) respectively. Here w r , u> y , and u> z are the vorticities in the x, y. 
and 2 directions, respectively, with the mean vorticity subtracted out. The angular 
brackets denote averaging over x-z planes. The agreement of the DM as well as 
the DLM(+ ) predictions with the filtered DNS is seen to be very good. When the 
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y/b 

FIGURE 4b. The mean cross-stream intensity of turbulence in self-similar coordi- 
nates using DLM(+) -f ; DM ; No model ; filtered DNS • 





// V 


FIGURE 4c. The mean spanwise intensity of turbulence in self-similar coordinates 
using DLM(+) +; DM ; No model ; filtered DNS • 



y/b 

FIGURE 4d. The mean turbulent stress uv in self-similar coordinates using 
DLM(+) +; DM ; No model ; filtered DNS • 
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IGURE 5a. The mean intensity of streamwise vorticity in self-similar coordinates 
sing DLM(+) +; DM ; No model x; filtered DNS • 



y/b 

FIGURE 5B . The mean cross-stream intensity of vorticity in self-similar coordinates 
using DLM(+) +; DM ; No model x; filtered DNS • 
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Figure 5c. The mean spanwise intensity of vorticity in self-similax coordinates 
using DLM(+) +; DM ; No model x; filtered DNS • 
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FIGURE 5d. The mean vorticity product uj x u) y in self-similar coordinates using 
DLM(+) +; DM ; No model x; filtered DNS • 
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model is turned off the agreement with the filtered DNS is seen to be very poor. The 
magnitudes of the enstrophy components are about four times the corresponding 
filtered DNS levels. Here one might ask if it is reasonable or useful to use an LES 
to predict vorticity statistics since it is known that small scales not resolved by LES 
Eire the primary contributors to enstrophy. Indeed, (a > 2 ), and (u> 2 z ) for the 

filtered DNS are about a fifth of their levels in the unfiltered DNS. However, good 
prediction of vorticity statistics is important because these statistics are a sensitive 
measure of the scales close to the threshold of resolution of the LES. The fact that 
even vorticity statistics are captured by the LES suggests that all of the resolved 
scales and not just the lowest wavenumber modes are faithfully represented in the 
simulation. Thus, we use vorticity statistics eis a “quality indicator” of the LES 
rather than as a quantity of practical importance to the user. 

In Figs. 4 and 5 it is apparent that the self-similar collapse is not perfect but 
that there is a systematic variation between the curves at different times in the 
simulation, even when scaled in self-similar variables. This is the case not only for 
the LES, but also for the filtered DNS. This is an artifact of the filtering procedure 
itself and can be understood in the following way. The flow evolves self-similarly at 
constant Reynolds number Ret, = b(6U m ) / v (see Fig. 2) in the self-similar region but 
the length scales increase in time. Thus, as the flow evolves, the energy spectrum 
shifts to the left without changing form. Since the grid size is held fixed, this 
implies that more and more of the energy becomes ‘resolved’ as the spectrum shifts 
to lower wavenumbers past k — 27r/A. Therefore the resolved part of the second- 
order statistics increases with time. This is precisely what is observed in the filtered 
DNS and LES data and is responsible for the systematic increasing trend during 
the self-similar period. 

In addition to obtaining quantitative predictions, one also hopes to gain some 
qualitative understanding of the large-scale flow structures from an LES. Thus, it 
is of interest to see if the model is able to generate structures that look realistic. 
As an example a typical contour plot of the u-velocity is presented in Fig. 6 over an 
x - y plane. It is seen that Fig. 6(C) (LES with model) bears an overall resemblance 
to Fig. 6(B) (filtered DNS) in the sense that it has a similar number of ‘eddies’ of 
approximately similar size and shape. However, Fig. 6(D) (LES without model) 
looks qualitatively different from Fig. 6(B) in the sense that it has a profusion of 
poorly resolved small-scale structures. A similar statement can be made about 
the other flow variables. The times at which the contours are shown in Fig. 6 
for the DNS and LES do not correspond exactly, but they are close, varying from 
t(6U 0 ) 2 /n « 62.4 to 66.3, and are in the developed region (see Fig. 1). 

It may appear that even though the no-model case (Fig. 6(D)) has far too many 
small-scale fluctuations compared to the filtered DNS (Fig. 6(B)), it does resemble 
somewhat the full DNS of Fig. 6(A). That this is not the case becomes clear on 
examination of the energy spectrum. Fig. 7 shows the one-dimensional spectrum 

q 2 (k z ) = y^(|ii(fc r , k z )\ 2 + |v(fci, k z )\ 2 + |u;(fc x , fc 2 )| 2 ) (H) 

k, 

at the plane y = 0 for the same fields whose u-contour plots are shown in Fig. 6. 



FIGURE 6. A typical iso- velocity contour plot for component v in the x 
for (A) DNS (B) Filtered DNS (C), DLM(+) and (D) No model. 
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Figure 7. The ‘(^-spectrum’ (see text) for LES using DLM(+ ) ! using no 

model ; filtered DNS • ; DNS . 

Here ‘ ~ ’ denotes Fourier-transform in the x — z plane. The filtered DNS is very 
close to, but slightly below (on account of filtering in the y-direction), the full DNS, 
up to the maximum k z represented in the LES. The LES matches the filtered DNS 
closely. However, for the no-model case, energy piles up at the high wavenumbers 
because of the lack of a dissipation mechanism, and this results in a ‘flat’ rather 
than decaying energy spectrum. The small-scale fluctuations seen in Fig. 6(D) 
are a manifestation of this unphysical ‘pile-up’ of energy and have no relation to 
the true fine structure seen in the highly resolved DNS of Fig. 6(A). Indeed, it is 
quite impossible to reproduce the fine structure of the DNS with the vastly reduced 
number of modes in an LES, and the ‘filtered’ DNS is the ideal limit one can hope 
to achieve. 

In summary, mean velocity profiles plotted in self-similar coordinates are very in- 
sensitive to the choice of subgrid models. The prediction of the self-similar growth of 
the wake width is improved by the subgrid model, but the results with no model are 
nevertheless tolerable. Second-order velocity and vorticity statistics are predicted 
very well by both the DM and DLM(+), but the predictions of these statistics 
without the model are very poor. The flow structures in the LES have a strong 
visual resemblance to those of the corresponding filtered DNS, but this is not the 
case if the LES is performed with no subgrid model. The LES represents a very 
significant saving in CPU time over the corresponding DNS. The results presented 
here suggest that LES can provide accurate predictions when information related 
to small-scale structures is not required. 
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3. Future plans 

In generating the results presented in this report, the “test filter” was imple- 
mented only in the x and 2 directions. This was done as a first step because filtering 
only in the homogeneous directions is easiest to implement in the code. However, 
there are some serious difficulties associated with this. In situations where the grid 
spacings in all three directions are not the same, there is no unique way of defining 
the “grid- filter width” A. Some possible choices are 

(A) A = (ArAyA*) 1 / 3 

(B) A = (A 2 + A 2 + A 2 ) 1 / 2 

(C) A=max [A I ,A !M A Z ] 

(D) A = (AjAj) 1 / 2 . 

The choice (D) may be thought of as a natural modification of (A) when A y <C 
A*, A z- In this case (A) would give an unreasonably small length scale. 

Now, if the filtering is done in all three directions, then A* = 2A X , A y = 2A # , and 
= 2A Z so that all the possible choices (A), (B), (C), and (D) give the same value 
for the filter- width ratio A/ A — 2. It is easily shown that only the combination C A? 
is computed in the dynamic model and that the grid and test filter-widths enter 
only as the ratio A/A. Thus when the model is properly implemented with a 3D- 
filter, it is unaffected by the choice of filter-width definition. In fact, any filter width 
A = f(A z ,A y ,A z ) where / is a homogeneous function (i.e. f(aA x ,bA y ,cA z ) = 
abcf(A x ,A y ,A z )) will yield the same filter-width ratio. This is no longer true if 
the filtering is only done in z — 2 planes as in the current simulations. In this case 
~ 2A X , Ay = Ay, and A z — 2A Z , so that (C) and (D) give A/A = 2. (A) 
gives A/A = 2 2 / 3 , and (B) gives a result that depends on the aspect ratio of the 
grid. In the simulations presented here we have chosen A/ A = 2, but the results 
change significantly if an alternate value for this ratio is used. These results should 
therefore be regarded as preliminary, and more careful tests using full 3D filtering 
need to be done before they can be considered reliable. 

We would like to test two other versions of the dynamic localization model viz. 
DLM(k) and DLM(S). The first one accounts for backscatter by means of a budget 
equation for the subgrid kinetic energy (Ghosal et al., 1992, 1994) while the second 
regards backscatter as a stochastic forcing. Apart from being able to represent 
backscatter (which may or may not be a significant effect), the DLM(k) has the 
additional advantage that it allows one to compute the full subgrid-stress tensor 
instead of simply the deviatoric part. This makes it possible to determine the 
resolved pressure, a quantity that cannot be determined if only the deviatoric part 
of the stress is known. 

Both DNS and experiments on plane wakes show a range of growth rates that 
seem to be sensitive to initial conditions (Moser and Rogers, 1994). It has been 
proposed that this could be due to the existence of non-unique self-similar states, 
any one of which can be selected in a given realization depending on the initial 
conditions (George, 1989). In order to investigate such possible dependence on 
initial conditions, Moser and Rogers (1994) amplified the 2D components of the 
initial velocity field in their simulation. It has been found that it is possible to 
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significantly alter the growth rate by such “2D forcing”. We would like to check if 
LES is able to predict correctly the growth rates in such forced wakes. If it does, 
then LES can be used as a research tool to test whether alternate self-similar states 
are indeed sustained. This requires long-time simulations that are prohibitively 
expensive using current DNS. 

We would like to thank Dr. Parviz Moin for his critical comments on an earlier 
version of this manuscript. 
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